cost <- read_csv("../data/MacroTrends_Data_Download.csv")
## Rows: 911 Columns: 3
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (1): date
## dbl (2): cost, nominal
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
# head(cost)
cost_ts <- cost %>% select(cost = cost) %>% ts(, start = c(1946, 1), frequency = 12)
ts_plot(cost_ts, 
        title = "World Yearly Oil Cost ",
        Ytitle = "Cost  in EJ",
        Xtitle = "Years",
        slider = TRUE)
dygraph(cost_ts, 
        main  = "World Yearly Oil Cost ",
        ylab  = "Cost  in EJ",
        xlab  = "Years") %>% 
  dyRangeSelector()
ts_decompose(cost_ts)
ts_cor(cost_ts)
ts_seasonal(cost_ts, type = "normal")
dy <- diff(diff(cost_ts))
autoplot(dy)

# autoplot(diff(cost_ts))
ggseasonplot(dy)

ggsubseriesplot(dy)

# Forecasting applications
# Setting training and testing partitions
cost_ts1 <- ts_split(ts.obj = cost_ts, sample.out = 12)
train <- cost_ts1$train
test <- cost_ts1$test

# Forecasting with auto.arima
library(forecast)
md <- auto.arima(train, d=1 )
summary(md)
## Series: train 
## ARIMA(1,1,0)(0,0,2)[12] 
## 
## Coefficients:
##          ar1     sma1     sma2
##       0.2000  -0.0715  -0.0605
## s.e.  0.0328   0.0341   0.0335
## 
## sigma^2 estimated as 21.85:  log likelihood=-2657.62
## AIC=5323.24   AICc=5323.29   BIC=5342.44
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.03426171 4.664044 2.533212 -0.1618761 4.375093 0.2565892
##                      ACF1
## Training set 0.0004869705
fc <- forecast(md, h = 12)
# time(fc)
# cycle(fc)
# class(fc)
# Plotting actual vs. fitted and forecasted
test_forecast(actual = cost_ts, forecast.obj = fc, test = test)
plot_forecast(fc)
# Run horse race between multiple models
methods <- list(ets1 = list(method = "ets",
                            method_arg = list(opt.crit = "lik"),
                            notes = "ETS model with opt.crit = lik"),
                ets2 = list(method = "ets",
                            method_arg = list(opt.crit = "amse"),
                            notes = "ETS model with opt.crit = amse"),
                arima1 = list(method = "arima",
                              method_arg = list(order = c(2,1,0)),
                              notes = "ARIMA(2,1,0)"),
                arima2 = list(method = "arima",
                              method_arg = list(order = c(2,1,2),
                                                seasonal = list(order = c(1,1,1))),
                              notes = "SARIMA(2,1,2)(1,1,1)"),
                hw = list(method = "HoltWinters",
                          method_arg = NULL,
                          notes = "HoltWinters Model"),
                tslm = list(method = "tslm",
                            method_arg = list(formula = input ~ trend + season),
                            notes = "tslm model with trend and seasonal components"))

# Training the models with backtesting
md <- train_model(input = cost_ts,
                  methods = methods,
                  train_method = list(partitions = 6, 
                                      sample.out = 12, 
                                      space = 3),
                  horizon = 12,
                  error = "MAPE")
## Warning in (function (x, alpha = NULL, beta = NULL, gamma = NULL, seasonal =
## c("additive", : optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## Warning in (function (x, alpha = NULL, beta = NULL, gamma = NULL, seasonal =
## c("additive", : optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
## # A tibble: 6 x 7
##   model_id model       notes avg_mape avg_rmse `avg_coverage_8~ `avg_coverage_9~
##   <chr>    <chr>       <chr>    <dbl>    <dbl>            <dbl>            <dbl>
## 1 arima1   arima       ARIM~    0.334     17.3            0.75             0.917
## 2 arima2   arima       SARI~    0.344     17.2            0.75             0.917
## 3 ets2     ets         ETS ~    0.359     18.4            0.5              0.75 
## 4 ets1     ets         ETS ~    0.365     18.1            0.528            0.806
## 5 hw       HoltWinters Holt~    0.370     18.7            0.542            0.889
## 6 tslm     tslm        tslm~    0.784     32.7            0.514            0.917
plot_model(md)
Cost_df <- data.frame(year = floor(time(cost_ts)), month = cycle(cost_ts),
   cost = as.numeric(cost_ts))
plot_forecast(fc)
# fc$model$model$# %>%  str()
forcast <- fortify(fc) 
cost_df <- data.frame(date  = forcast$Index, cost = c(forcast [ 1:899, 3], forcast [ 900:nrow(forcast), 4]))
cost_df %>% ggplot(aes(x = date, y = cost)) +
  geom_line()